
#### Produces Figure A4 ######

install.packages("MCMCpack")

library(MCMCpack)


#setwd("~/Dropbox/Job Market Paper/IO Revisions/Replication_Files/Replication_Appendix/A_Fig_4")



data_q<-read.csv(file="data_fig_A4.csv")


num<-tapply(data_q$ID,data_q$Year,length)


years<-seq(1105,1790,by=5)



############# This Takes a Long Time To Run ######

output2<-list()

priorsa<-c(10,1.5,1/2)
priorsbM<-c(1,1,1)

for(i in 1:3){
  
  out1<-MCMCpoissonChange(num~1,m=1,c0=1,d0=1,a=139/20 * i,b=1,mcmc=100000,burnin=100000,marginal.likelihood = "Chib95")
  
  out2<-MCMCpoissonChange(num~1,m=2,c0=1,d0=1,a=139/30 * i,b=1,mcmc=100000,burnin=100000,marginal.likelihood = "Chib95")
  
  out3<-MCMCpoissonChange(num~1,m=3,c0=1,d0=1,a=139/40 * i,b=1,mcmc=100000,burnin=100000,marginal.likelihood = "Chib95")
  
  out4<-MCMCpoissonChange(num~1,m=4,c0=1,d0=1,a=139/50 * i,b=1,mcmc=100000,burnin=100000,marginal.likelihood = "Chib95")
  
  out5<-MCMCpoissonChange(num~1,m=5,c0=1,d0=1,a=139/60 * i,b=1,mcmc=100000,burnin=100000,marginal.likelihood = "Chib95")
  
  out6<-MCMCpoissonChange(num~1,m=6,c0=1,d0=1,a=139/70 * i,b=1,mcmc=100000,burnin=100000,marginal.likelihood = "Chib95")
  
  out7<-MCMCpoissonChange(num~1,m=7,c0=1,d0=1,a=139/80 * i,b=1,mcmc=100000,burnin=100000,marginal.likelihood = "Chib95")
  
  BF<-BayesFactor(out1,out2,out3,out4,out5,out6,out7)
  
  output2[[i]]<-list(out1,out2,out3,out4,out5,out6,out7,BF)
  
  
}  


print("This is the Bayes Factor Selection - Chooses Model 1")

print(BF)


trans<-list()

p<-runif(10,0,1)



for(i in 1:10){
  
trans[[i]]<-matrix(NA,2,2)  
  
trans[[i]][1,1]<-p[i]
trans[[i]][1,2]<-1-p[i]
trans[[i]][2,1]<-0
trans[[i]][2,2]<-1
  
}



Bs<-rnorm(10,log(mean(data_q$logArea)),sd(data_q$logArea))

   
   
   
   
probs<-attributes(out1)[13]



pdf(file="Figure_AF.pdf")

plot(years,num[2:139],type="l",ylab="Number of States", xlab="Time")

par(new=T)

plot(years,probs[[1]][2:139,1],lty=2,axes=FALSE,ylim=c(0,1),type="l",ylab="",xlab="")

par(new=T)

plot(years,probs[[1]][2:139,2],lty=3,axes=FALSE,ylim=c(0,1),type="l",ylab="",xlab="")

axis(4,at=c(0,1),cex=.5)
mtext("P(S = k| Y)", 4, line=1)




p<-runif(10,0,1)



for(i in 1:10){
  
  trans[[i]]<-matrix(NA,2,2)  
  
  trans[[i]][1,1]<-p[i]
  trans[[i]][1,2]<-1-p[i]
  trans[[i]][2,1]<-0
  trans[[i]][2,2]<-1
  
}


dev.off()


save.image(file="Fig_A4.Rdata")


